A  HOUGH  TRANSFORM  GLOBAL  PROBABILISTIC  APPROACH 
TO  MULTIPLE  SUBJECT  DIFFUSION  MRI  TRACTOGRAPHY 


By 

Iman  Aganj,  Christophe  Lenglet 
Neda  Jahanshad,  Essa  Yacoub 
Noam  Harel,  Paul  M.  Thompson 

and 

Guillermo  Sapiro 


IMA  Preprint  Series  #  2305 

(April  2010) 


INSTITUTE  FOR  MATHEMATICS  AND  ITS  APPLICATIONS 

UNIVERSITY  OF  MINNESOTA 

400  Lind  Hall 
207  Church  Street  S.E. 

Minneapolis,  Minnesota  55455-0436 

Phone:  612-624-6066  Fax:  612-626-7370 
URL:  http://www.ima.umn.edu 


Form  Approved 
OMB  No.  0704-0188 


Report  Documentation  Page 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

APR  2010 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

4.  TITLE  AND  SUBTITLE 

A  Hough  Transform  Global  Probabilistic  Approach  to  Multiple-Subject 
Diffusion  MRI  Tractography 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Minnesota, Institute  for  Mathematics  and  Its 

Application, 207  Church  Street  SE, Minneapolis, MN, 55455-0436 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

A  global  probabilistic  fiber  tracking  approach  based  on  the  voting  process  provided  by  the  Hough 
transform  is  introduced  in  this  work.  The  proposed  framework  tests  candidate  3D  curves  in  the  volume 
assigning  to  each  one  a  score  computed  from  the  diffusion  images,  and  then  selects  the  curves  with  the 
highest  scores  as  the  potential  anatomical  connections.  The  algorithm  avoids  local  minima  by  performing 
an  exhaustive  search  at  the  desired  resolution.  The  technique  is  easily  extended  to  multiple  subjects 
considering  a  single  representative  volume  where  the  registered  high-angular  resolution  diffusion  images 
(HARDI)  from  all  the  subjects  are  non-linearly  combined,  thereby  obtaining  population-representative 
tracts.  The  tractography  algorithm  is  run  only  once  for  the  multiple  subjects,  and  no  tract  alignment  is 
necessary.  We  present  experimental  results  on  HARDI  volumes,  ranging  from  a  1.5T  physical  phantom  to 
7T  and  4T  human  brain  and  7T  monkey  brain  datasets. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

25 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


A  Hough  Transform  Global  Probabilistic  Approach  to  Multiple- 
Subject  Diffusion  MRI  Tractography 


Iman  Aganj,a  Christophe  Lenglet,ab  Neda  Jahanshad,c  Essa  Yacoub,b  Noam  Harel,b  Paul  M.  Thompson,0 

and  Guillermo  Sapiroa 


a  Department  of  Electrical  and  Computer  Engineering,  University  of  Minnesota,  Minneapolis,  200  Union 
St.  SE,  MN  55455,  USA 

b  Center  for  Magnetic  Resonance  Research,  University  of  Minnesota,  2021  Sixth  Street  SE,  Minneapolis, 
MN  55455,  USA 

c  Laboratory  of  Neuro  Imaging,  University  of  California-Los  Angeles,  School  of  Medicine,  635  Charles 
Young  Drive  South,  Suite  225,  Los  Angeles,  CA  90095-7334,  USA 


Corresponding  author: 

Iman  Aganj 

200  Union  St.  SE,  Minneapolis,  MN  55455,  USA 
Tel:  +1  (612)  202-3019 

Email:  iman@umn.edu 


1 


ABSTRACT 


A  global  probabilistic  fiber  tracking  approach  based  on  the  voting  process  provided  by  the  Hough 
transform  is  introduced  in  this  work.  The  proposed  framework  tests  candidate  3D  curves  in  the  volume, 
assigning  to  each  one  a  score  computed  from  the  diffusion  images,  and  then  selects  the  curves  with  the 
highest  scores  as  the  potential  anatomical  connections.  The  algorithm  avoids  local  minima  by  performing 
an  exhaustive  search  at  the  desired  resolution.  The  technique  is  easily  extended  to  multiple  subjects, 
considering  a  single  representative  volume  where  the  registered  high-angular  resolution  diffusion  images 
(HARDI)  from  all  the  subjects  are  non-linearly  combined,  thereby  obtaining  population-representative 
tracts.  The  tractography  algorithm  is  run  only  once  for  the  multiple  subjects,  and  no  tract  alignment  is 
necessary.  We  present  experimental  results  on  HARDI  volumes,  ranging  from  a  1.5T  physical  phantom  to 
7T  and  4T  human  brain  and  7T  monkey  brain  datasets. 


Keywords:  Tractography,  diffusion-weighted  magnetic  resonance  imaging  (DWI),  Hough  transform, 

orientation  distribution  function  (ODF),  population  studies. 
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1.  Introduction 


Understanding  the  connectivity  between  different  areas  of  the  brain  is  essential  in  studying  brain 
function  and  development.  Diffusion-weighted  magnetic  resonance  imaging  (DWI)  provides,  through 
tractography,  a  unique  in-vivo  quantitative  measurement  of  the  brain’s  anatomical  connectivity.  In 
addition  to  its  benefits  in  neurosurgical  planning,  DWI  tractography  has  considerable  clinical  importance 
by  noninvasively  quantifying  changes  in  the  white  matter  connectivity  at  different  stages  of  diseases  or 
development.  Moreover,  it  can  be  used  to  segment  fiber  bundles  of  the  central  nervous  system,  or  in  tract- 
based  statistical  analysis  of  scalars  such  as  the  fractional  anisotropy  (FA).  Performing  tractography  in 
multiple  subjects  is  invaluable  for  population  studies  and  creating  fiber  bundle  atlases. 

DWI  provides  local  information  about  the  fiber  orientation  by  measuring  the  tissue  water  diffusion,  in 
vivo ,  assuming  a  high  correlation  between  the  fiber  and  diffusion  orientations.  However,  there  is  no 
unique  solution  as  to  how  to  integrate  these  voxel-scale  local  orientations  to  infer  global  connectivity. 
Early  fiber  tractography  algorithms,  known  as  streamline  methods,  are  based  on  following  the  principal 
diffusion  orientation  (Basser  et  al.,  2000;  Conturo  et  al.,  1999;  Jones  et  al.,  1999;  Lazar  et  al.,  2003;  Mori 
et  al.,  1999).  Despite  their  simplicity,  these  methods  are  prone  to  cumulative  errors  caused  by  noise, 
partial  volume  effects,  and  discrete  integration,  and  have  difficulty  in  distinguishing  fiber  crossing  and 
kissing  mostly  due  to  the  fact  that  the  entire  diffusion  information  is  not  globally  used  and  integrated. 
This  led  to  the  development  of  other  successful  approaches,  including  probabilistic  techniques  (Behrens  et 
al.,  2007;  Bjornemo  et  al.,  2002;  Descoteaux  et  al.,  2009;  Friman  et  al.,  2006;  Jones,  2008;  Lazar  & 
Alexander,  2005;  Parker  et  al.,  2003), global  techniques  based  on  front  propagation  (Campbell  et  al.,  2005; 
Jackowski  et  al.,  2005;  Parker  et  al.,  2002;  Pichon  et  al.,  2005;  Prados  et  al.,  2006;  Tournier  et  al.,  2003), 
simulation  of  the  diffusion  process  or  fluid  flow  (Batchelor  et  al.,  2001;  Hageman  et  al.,  2009;  Hagmann 
et  al.,  2003;  Kang  et  al.,  2005;  ODonnell  et  al.,  2002;  Yortik  et  al.,  2005),  DWI  geodesic  computations 
(Jbabdi  et  al.,  2008;  Lenglet  et  al.,  2009a;  Melonakos  et  al.,  2007;  Pechaud  et  al.,  2009),  graph  theoretical 
techniques  (Iturria-Medina  et  al.,  2007;  Sotiropoulos  et  al.,  2009;  Zalesky,  2008),  spin  glass  models 
(Fillard  et  al.,  2009;  Mangin  et  al.,  2002),  and  Gibbs  tracking  (Kreher  et  al.,  2008).  Generally  speaking, 
for  virtually  every  tractography  method,  a  particular  putative  subset  of  all  possible  curves  is  implicitly 
considered  from  which  the  resulting  tracts  are  chosen  according  to  some  criteria.  This  subset,  although 
different  depending  on  the  particular  sweeping  strategy,  virtually  never  coincides  with  the  universal  set  of 
all  possible  curves  representing  the  candidate  tracts.  The  closer  the  subset  is  to  the  universal  set  of  curves, 
the  more  accurate  we  expect  the  results  to  be.  For  a  recent  thorough  discussion  on  different  tractography 
techniques,  see  (Behrens  &  Jbabdi,  2009). 
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Prior  approaches  for  multi-subject  tractography  are  typically  based  on  the  post  processing  of 
tractography  results  from  individual  subjects  (El  Kouby  et  al.,  2005;  Jbabdi  et  al.,  2009;  Leemans  et  al., 
2006;  Maddah  et  al.,  2006;  O'Donnell  &  Westin,  2007;  Voineskos  et  al.,  2009;  Wakana  et  al.,  2004). 
These  methods  generally  require  aligning  the  tracts  and  mapping  them  into  a  common  fiber  coordinate 
system,  which  is  challenging  due  to  the  large  number  of  high-dimensional  fiber  trajectories  per  subject 
and  the  lack  of  clearly  defined  criteria  for  aligning  curves  and  particularly  tracts. 

In  this  work,  we  present  a  global  probabilistic  approach  inspired  by  the  voting  procedure  provided  by 
the  popular  Hough  transform  (Duda  &  Hart,  1972;  Gonzalez  &  Woods,  2008).  Our  proposed  tractography 
algorithm  essentially  tests  candidate  3D  curves  in  the  volume,  assigning  a  score  to  each  of  them,  and  then 
returning  the  curves  with  the  highest  scores  as  the  potential  anatomical  connections.  The  score  is 
accordingly  derived  from  the  DWI  data.  Being  an  exhaustive  search,  this  proposed  algorithm  avoids 
entrapment  in  local  minima  within  the  discretization  resolution  of  the  parameter  space.1  Furthermore,  the 
specific  definition  of  the  candidate  tract  score  has  the  desired  effect  of  attenuating  the  noise  through  the 
integration  of  the  real-valued  local  votes  derived  from  the  diffusion  data.  We  also  introduce  a 
simultaneous  multi-subject  tractography  technique  which  takes  as  input  a  single  representative  volume  - 
where  the  HARDI  data  from  all  the  (registered)  subjects  are  non-linearly  integrated  -  and  generates 
population-representative  tracts.  The  multi-subject  tractography  algorithm  is  run  only  once,  and  no  tract 
alignment  is  necessary.  We  present  experimental  results  on  HARDI  volumes  such  as  a  biological 
phantom  dataset  acquired  at  1.5T,  a  monkey  brain  dataset  acquired  at  7T,  and  a  number  of  human  brain 
datasets  acquired  at  4T  and  7T.2 

In  Sec.  2  we  present  the  proposed  algorithm  in  detail.  Experimental  results  are  presented  in  Sec.  3, 
and  Sec.  4  concludes  with  a  review  of  the  contributions.  Additional  implementation  details  are  provided 
in  the  Appendix. 


2.  Methods 

We  first  randomly  generate  a  sufficiently  high  number  of  initial  seed  points  inside  a  brain  mask  or  a 
region  of  interest.  From  each  initial  point,  we  consider  as  many  passing  curves  as  desired,  based  on  the 
expected  resolution  and  available  computational  resources  (Fig.  1,  left).  A  score  is  computed  for  each 

1  Please  note  that  many  prior  approaches  may  as  well  be  modified  to  perform  efficient  exhaustive  searches. 

2  This  paper  extends  our  previous  conference  versions  for  single  and  multiple  subjects  (Aganj  et  al.,  2009a, b).  In 

particular,  we  provide  more  implementation  details  and  additional  validation  and  comparisons. 
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curve,  and  the  one  (or  ones)  with  the  maximum  score  is  then  chosen  as  the  best  curve(s)  representing  the 
fiber  bundle  passing  through  that  seed  point  (Fig.  1,  middle  &  right).3  This  process  is  detailed  in  the 
following  subsections. 


2.7.  Curve  parameterization 

We  parameterize  the  3D  curves  by  the  arc  length  s ,  with  5  =  0  corresponding  to  the  seed  point  (Fig. 
2).  The  unit  tangent  vector  of  the  curve  is  identified  at  each  point  by  standard  polar  coordinates  0(s )  and 
00): 

(sin0(s)  cos  0(s)\ 

sin0(s)  sin0(s)  j.  (1) 

cos  0(s)  ) 

In  our  proposed  model,  we  consider  simple  polynomial  approximations  of  these  two  angles  with 
respect  to  the  arc  length: 

1 v 

00)  =  ^  CLksk, 

k= 0 

N 

000  =  ^  bksk, 

k= 0 

where  N  is  the  polynomial  order  (different  orders  for  6  and  0  can  be  considered  if  desired).  The  Hough- 
inspired  process  will  be  used  to  select  the  best  possible  coefficients  ak  and  bk  (and  the  two  additional 
parameters  introduced  below),  based  on  the  available  diffusion  data.  In  addition,  two  extra  parameters  L_ 
and  L+  determine  the  partial  (Euclidean)  lengths  of  the  curve  on  each  side  of  the  seed  point  (Fig.  2),  with 
0  <  L_  +  <  Lmax  where  the  constant  Lmax  is  chosen  as  the  maximum  expected  curve  (fiber)  length.  Each 
curve  initiated  from  the  seed  point  x0  is  then  represented  using  d  =  2N  +  4  unique  parameters, 
{ a0 , ... ,  aN,  b0, ... ,  bN,  L_,  L+},  and  explicitly  computed  by  integrating  the  tangent  vector: 


(2) 

(3) 


3  As  customary  in  probabilistic  techniques,  several  candidate  curves  may  as  well  be  selected  per  seed  point,  each 
carrying  a  score. 
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(4) 
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x(s)  =  x0  + 


s  E  [—L_,L+\. 


2.2.  Fiber  score  computation 

A  score,  intended  to  estimate  the  log-probability  of  the  existence  of  a  fiber,  is  assigned  to  each 
possible  3D  curve  passing  through  a  seed  point  x0-  In  this  work,  the  score  is  defined  as 


T 

Sx0(a0,  ...,aN,b0,  ...,bN,L_,L+)  ~  J  (log P(x(s),t(s))  +  A)ds 


-L_ 


(5) 


—  log  P(x(s),  t(s))  ds  +  A(L_  +  L+). 

j 


-L_ 


The  expression  P(x,  f)dfl  represents  the  probability  for  the  point  x  to  be  located  inside  a  fiber  bundle 
passing  in  the  direction  £  through  the  infinitesimal  solid  angle  dfl.4  The  constant  A,  which  is  used  to 
compensate  for  the  absence  of  dfl  in  the  integral,  can  also  be  interpreted  as  a  prior  on  the  length  of  the 
fiber  bundles,  as  choosing  a  larger  A  favors  the  score  of  longer  curves.5 

P(x,  £)  can  be  computed  using  the  conditional  probability  formula  as, 

P(x,  £)  =  P(x)P(£|x).  (6) 

The  prior  probability  of  the  existence  of  a  fiber  at  the  point  x,  P(x),  is  considered  to  be  equal  to  either 
the  fractional  anisotropy  (FA)  or  generalized  fractional  anisotropy  (GFA)  inside  the  brain  tissue,  and  zero 
outside  the  brain  mask  and  inside  the  cerebrospinal  fluid.  This  comes  from  the  assumption  that  the  more 
anisotropic  a  region  is,  the  more  likely  a  fiber  bundle  may  be  passing  through  that  region.  In  addition,  as 
long  as  no  further  constraints  or  selections  are  provided  by  the  user,  the  initial  seed  points  are  chosen 


4  The  curve  score  is  defined  by  assuming  that  adjacent  voxels  are  independent,  making  it  only  an  approximation  of 
the  true  fiber  log-probability.  The  incorporation  of  spatial  coherence  and  continuity  is  the  subject  of  future  research. 

5  Not  including  A  may  cause  the  zero-length  curve  (L_  =  L+  =  0)  to  gain  the  maximum  score,  due  to  the  potentially 

negative  values  of  the  logarithm. 
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randomly  with  a  spatial  probability  distribution  proportional  to  P(x).  Other  choices  for  P(x),  such  as  the 
white  matter  complexity  introduced  in  (Haro  et  al.,  2008),  are  also  possible. 

Next,  assuming  that  a  fiber  is  actually  passing  through  the  point  x,  the  probability  for  it  to  be  in  the 
direction  f,  i.e.  P(t |x),  is  derived  from  the  orientation  distribution  function  (ODF)  at  each  voxel  in  the 
volume.  Computed  from  various  DWI  modalities,  the  ODF  is  defined  as, 

n  00 

0DF(u )  :=  PDF(ru)r2dr,  (7) 

Jo 

which  is  the  integration  of  the  PDF(r )  (the  spatial  probability  density  function  of  the  diffusion  of  water 
after  a  certain  amount  of  time),  in  a  cone  of  constant  solid  angle  (CSA)  in  the  direction  of  the  unit  vector 
u.  In  Diffusion  Spectrum  Imaging  (DSI)  (Wedeen  et  al.,  2005),  PDF(r )  is  available  on  a  discrete 
Cartesian  grid,  and  therefore  the  ODF  is  directly  computed  from  the  above  formula.  In  the  case  of  the 
Diffusion  Tensor  Imaging  (DTI)  (Basser  et  al.,  1994),  the  ODF  is  computed  by  integrating  the  3D  normal 
distribution, 


ODF  (u)  =  f - \ - -e_l™T£>_1™r2dr  = - j-2 - -  (8) 

mDTI  o  (2tt)2|D|2  4tt|D|2(utD-1u)2 

where  D  is  proportional  to  the  estimated  diffusion  tensor. 

In  this  work,  however,  we  use  Q-ball  Imaging  (QBI)  (Tuch,  2004),  which  is  a  popular  HARDI 
reconstruction  method  proven  successful  in  resolving  multiple  intravoxel  fiber  orientations.  The  original 
ODF  expression  in  QBI  does  not  include  the  Jacobian  factor  r2,  creating  the  need  for  post-processing 
such  as  artificial  sharpening.  Here  we  use  the  normalized  and  dimensionless  ODF  estimator  in  QBI, 
derived  in  (Aganj  et  al.,  2010),  which  by  considering  the  factor  r2  computes  the  CSA-ODF, 

ODF  (u)  «T  +  ^_FRT(v2in(_in^M)]  (9) 

in  CSA-QBP  2  4?T  16tT2  [  b  \  S0  )) 

with  S(u)  and  S0  the  diffusion  signal  and  the  baseline  image  respectively,  and  FRT  and  the  Funk- 
Radon  transform  (Funk,  1916)  and  the  Laplace-Beltrami  operator  respectively.  This  ODF  reconstruction 
scheme  has  been  shown  to  outperform  the  original  QBI  by  improving  the  resolution  of  multiple  fiber 
orientations  (Aganj  et  al.,  2010),  and  producing  more  stable  and  consistent  GFA  (Fritzsche  et  al.,  2010). 
To  allow  sampling  in  any  desired  direction  f,  the  ODFs  were  approximated  in  the  real  and  symmetric 
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modified  spherical  harmonic  basis,  following  the  method  proposed  by  Descoteaux  et  al.  (2007)  for  the 
original  QBI,  and  subsequently  adapted  in  (Aganj  et  al.,  2010)  for  the  CSA-QBI. 

Putting  all  this  together,  and  using  for  instance  the  FA  as  P(x),  the  score  in  Eq.  (5)  thus  becomes, 

Sx0 (,a0'  —  >  aN> b0, bN, L_, L+)  :=  J  (log[ODFjCs)({(s))FA(x(s))]  +  X)ds,  (10) 

— L_ 

where  0DFf(5)(t(s))  stands  for  the  ODF  at  the  3D  position  x(s)  evaluated  in  the  direction  t(s),  with 
x(s)  and  t(s)  represented  via  the  polynomials  as  specified  in  equations  (1-4).  The  score  integral  in 
equations  (5)  and  (10)  has  the  additional  nice  effect  of  attenuating  the  additive  noise  in  the  data  through 
summation. 


2.3.  Hough  transform 

As  discussed  in  Sec.  2.1,  every  curve  starting  from  a  particular  seed  point  is  presented  as  a  point  in  a 
d-dimensional  space,  with  d  —  2N  +  4  being  the  number  of  necessary  parameters.  In  theory,  we  would 
like  to  find  all  possible  curves  which  pass  through  the  seed  point  while  computing  their  scores,  to 
eventually  choose  the  one(s)  with  the  highest  score(s)  as  the  potential  fiber  tract(s)  passing  through  the 
seed  point.  However,  we  can  only  perform  such  an  exhaustive  search  within  a  finite  resolution,  by 
discretizing  Rd  and  assigning  discrete  values  to  the  curve  parameters  within  some  predefined  limits.6  The 
resulting  d-dimensional  array  of  curve  scores  is  often  called  the  Hough  transform  (Duda  &  Hart,  1972; 
Gonzalez  &  Woods,  2008)  of  the  data  with  respect  to  the  curves  passing  through  the  chosen  seed  point. 
This  can  be  seen  as  a  voting  process  where  the  voxels  cast  real-valued  votes  for  the  curves.  The  overall 
vote  is  the  integrand  of  the  score  integral  (Eq.  (5)  or  (10))  if  a  curve  passes  through  a  voxel,  and  zero 
otherwise. 

The  proposed  method  avoids  entrapment  in  local  minima  by  performing  an  exhaustive  search  in  the 
(discretized)  high-dimensional  space  of  the  curves.  Nevertheless,  the  discretization  resolution  of  the 
parameter  space  causes  the  algorithm  to  obtain  an  approximation  of  the  true  global  optimum  (which  is 
improved  by  increasing  the  resolution  as  desired).  To  alleviate  this  issue,  we  choose  the  best  curve  in  a 
multi-resolution  approach:  once  the  point  (in  the  parameter  space  Rd)  corresponding  to  the  curve  with  the 
highest  score  is  found  in  one  resolution  level,  the  neighborhood  of  that  point  is  discretized  again  with  a 

6  This  is  the  standard  procedure  in  the  Hough  transform,  where  the  accumulator  is  discretized. 
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higher  resolution  and  the  search  is  continued  at  the  next  level.  We  have  performed  our  experiments  using 
three  levels  of  resolution. 

This  concludes  the  description  of  the  proposed  technique  for  a  single  dataset.  We  now  show  how  this 
can  be  efficiently  extended  to  multiple  datasets  or  subjects. 


2.4.  Extension  to  multiple  subjects 

Here  we  extend  our  Hough  transform-based  global  approach  to  obtain  average  representative  tracts 
from  multiple  subjects  or  datasets.  We  perform  this  by  first  registering  the  HARDI  volumes,  using  either 
linear  transformation  or  more  sophisticated  algorithms  such  as  (Chiang  et  al.,  2008),  and  then  running  the 
algorithm  described  above  on  a  single  equivalent  volume  composed  of  the  voxel-wise  mean  ODF  and 
mean  FA  across  all  the  subjects.7  We  may  use  either  the  arithmetic  or  the  geometric  mean,  however,  the 
linearity  of  the  curve  score  (Eq.  (10))  with  respect  to  the  logarithms  of  the  ODF  and  FA  makes  the  use  of 
the  geometric  mean  more  appealing  (since  the  arithmetic  mean  of  the  logarithms  of  the  ODF  and  FA 
values  equals  the  logarithm  of  their  geometric  mean).  Hence,  we  reconstruct  the  effective  ODF  and  FA 
for  each  voxel  by  computing  the  geometric  mean  of  their  values  across  the  subjects, 

l 

(M 

UODF& 

i= 1 

(M  \  M 

PjFi4‘(x)J  ,  (12) 

where  the  superscripts  i  and  eq  indicate  respectively  the  ith  subject  (out  of  a  total  of  M)  and  the  equivalent 
subject.  We  eventually  use  the  equivalent  ODF  and  FA  volumes  in  the  single-subject  tractography 
algorithm,  thereby  running  it  only  once  for  all  the  subjects  and  avoiding  complications  due  to  curve 
(tracts)  registration. 


ODF^q(i) 


FAe%x) 


7 


This  could  be  interpreted  as  multiple  votes  per  voxel,  cast  by  each  corresponding  voxel  in  each  volume. 
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3.  Experimental  results  and  discussion 

3.1.  Results  for  single  subjects 

We  tested  our  method  on  various  HARDI  datasets,  for  each  of  which  the  FA  and  the  CSA-ODFs  were 
computed  as  explained  in  Sec.  2.2. 

We  first  used  the  biological  phantom  in  (Campbell  et  al.,  2005),  constructed  from  excised  rat  spinal 
cords  and  designed  to  have  crossing  tracts  (90  diffusion  images  at  b  =  1300  s/mm2,  1.5T).  We  computed 
the  tracts  from  200  seed  points,  using  three  different  bias  values  of  2  =  2.0,  2.5,  and  3.0  (see  Sec.  2.2),  and 
polynomials  of  order  N  =  3,  resulting  in  a  total  number  of  d  —  10  parameters  to  represent  the  candidate 
3D  curves  initiated  from  each  seed  point.  Figure  3  (top,  left)  shows  the  ODFs  superimposed  on  the  FA 
map,  and  the  rest  of  the  subfigures  show  the  tractography  results  using  different  values  for  X.  Increasing  X 
results  in  longer  curves  being  selected.  The  color  and  the  opacity  of  each  tract  (in  all  the  figures)  increase 
with  the  score,  from  transparent  blue  to  opaque  red. 

We  performed  additional  experiments  on  a  human  brain  HARDI  dataset  acquired  at  7T.  A  single 
refocused  2D  single  shot  spin  echo  EPI  sequence  was  used.  Image  parameters  were:  FOV:  192x192  mm2 
(matrix:  196x96)  to  yield  a  spatial  resolution  of  2x2x2  mm3,  TR/TE  4800/57  msec.,  acceleration  factor 
(GRAPPA)  of  2  and  6/8  partial  Fourier  were  used  along  the  phase  encode  direction.  Diffusion-weighted 
images  were  acquired  at  b  =  3000  s/mm2  with  256  directions,  along  with  31  baseline  images.  EPI  echo 
spacing  was  0.57  msec,  with  a  bandwidth  of  2895  Hz/Pixel.  Tracts  were  computed  from  1500  seed  points 
in  two  experiments,  using  polynomials  of  orders  N  —  2  and  N  —  3  (Fig.  4).  High  scoring  curves  are 
concentrated  in  major  fiber  bundles  such  as  corona  radiata,  corpus  callosum,  cingulum,  superior 
longitudinal  fasciculus,  and  arcuate  fasciculus.  Higher  polynomial  order  brings  more  flexibility  to  the 
curves,  resulting  in  them  being  spread  out  more  continuously  in  the  white  matter  regions  (e.g.  in  corona 
radiata).  A  3D  stereoscopic  rendering  of  the  results  is  shown  in  Fig.  5. 

Next,  we  used  the  monkey  brain  HARDI  dataset  introduced  in  (Lenglet  et  al.,  2009b)  to  test  the 
performance  of  our  method  on  specific  fiber  bundles.  An  anesthetized  Macaca  mulatta  monkey  was 
scanned  using  a  7T  MR  scanner  (Siemens)  equipped  with  a  head  gradient  coil  (80mT/m  G-maximum, 
200mT/m/ms)  with  a  diffusion  weighted  spin-echo  EPI  sequence.  Diffusion  images  were  acquired  at  b  = 
3000  s/mm2  (twice  during  the  same  session,  and  then  averaged)  over  100  directions  uniformly  distributed 
on  the  sphere.  We  used  TR/TE  of  4600/65  ms,  and  a  voxel  size  of  lxlxl  mm3.  We  computed  the  tracts 
from  1350  seed  points  uniformly  distributed  in  a  mask  containing  the  intersection  of  the  forceps  minor 
and  the  inferior  longitudinal  fasciculus,  using  the  polynomials  of  order  N  —  3.  Results  are  depicted  in 
Fig.  6.  A  fiber  density  map  was  created  by  counting,  at  each  voxel,  the  number  of  intersecting  curves 
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while  taking  into  account  their  respective  score.  A  three-dimensional  isosurface  was  then  generated  by 
thresholding  this  map  to  keep  the  most  relevant  connections.  It  is  as  well  presented  in  Fig.  6,  overlaid  on  a 
structural  MRI.  Major  tracts  including  the  splenium  fibers,  posterior  corona  radiata,  tapetum,  as  well  as 
the  inferior  fronto-occipital  and  longitudinal  fasciculi  -  including  the  optic  radiations  -  are  clearly 
identified.  Moreover,  fibers  of  the  optic  tract  are  recovered  until  they  reach  the  optic  chiasm. 


3.2.  Results  for  multiple  subjects 

We  used  our  multi-subject  tractography  algorithm  to  compute  mean  tracts  from  five  HARDI  datasets, 
introduced  in  (de  Zubicaray  et  al.,  2008),  each  acquired  from  a  different  healthy  young  adult.  Images  were 
acquired  using  a  4T  Bruker  Medspec  MRI  scanner.  Diffusion-weighted  images  were  acquired  using 
single-shot  echo  planar  imaging  with  a  twice-refocused  spin  echo  sequence  to  reduce  eddy-current 
induced  distortions.  Imaging  parameters  were:  23  cm  FOV,  TR/TE  6090/91.7  ms,  with  a  128x100 
acquisition  matrix.  Each  3D  volume  consisted  of  55  2-mm  thick  axial  slices  with  a  1.8x1. 8  mm2  in-plane 
resolution.  105  images  were  acquired:  11  with  no  diffusion  sensitization  (i.e.,  T2-weighted  b0  images)  and 
94  diffusion- weighted  images  (b  =  1159  s/mm2)  with  gradient  directions  evenly  distributed  on  the 
hemisphere.  Scan  time  was  approximately  14  minutes.  Images  were  corrected  for  motion  and  eddy 
current  distortions.  Each  subject’s  average  b0  image  was  aligned  to  a  group-specific  minimal  deformation 
template  (MDT)  using  a  nine-parameter  affine  transformation.  This  transformation  was  then  applied  to 
each  individual  DWI  and  gradient  directions  were  corrected  accordingly  for  ODF  calculations. 

We  combined  the  individual  datasets  into  two  equivalent  volumes,  using  the  geometric  and  the 
arithmetic  means  (see  Sec.  2.4).  We  tested  the  tractography  algorithm  on  both  equivalent  volumes,  and 
also  for  comparison,  on  two  of  the  five  individual  datasets.  In  each  experiment,  polynomials  of  order 
N  =  3  were  used  to  represent  3D  tracts  initiated  from  1500  seed  points.  Figure  7  (two  top  rows)  show  the 
mean  tracts  from  the  five  subjects  using  respectively  the  geometric  and  arithmetic  means.  The  two  bottom 
subfigures  show  tracts  from  individual  subjects.  As  Fig.  7  demonstrates,  combining  the  volumes  improves 
the  results  by  decreasing  the  number  of  false  positives  produced  by  noise.  In  addition,  the  fibers  are  less 
scattered  and  better  concentrated  in  major  fiber  bundles.  Note  particularly  how  corticospinal  tracts  are 
enhanced. 
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4.  Conclusion 


We  have  introduced  a  global  approach  for  single-  and  multi-subject  probabilistic  tractography,  based 
on  the  voting  process  provided  by  the  Hough  transform.  We  presented  experimental  results  on  a  physical 
phantom  and  brain  HARDI  datasets,  and  showed  that  using  this  approach,  data  from  multiple  subjects  can 
be  non-linearly  combined  and  exploited  to  obtain  population  statistics  and  more  accurate  tractography 
results.  The  incorporation  of  spatial  coherence  and  continuity  in  curve  score  computation  is  the  subject  of 
future  research. 
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Appendix 

In  this  appendix,  we  provide  additional  details  on  the  implementation  of  the  proposed  tractography 
technique. 

The  Hough  transform  is  often  used  in  global  optimization  problems  to  avoid  local  optimum  solutions. 
This  approach  is  however  characterized  by  its  high  computational  complexity,  given  that  in  a 
straightforward  implementation,  all  possible  solutions  must  be  tested  in  order  to  reconstruct  the  table  of 
scores.  This  issue  could  yet  be  alleviated  by  parallelizing  the  exhaustive  searches  at  the  seed  points,  as 
they  can  be  computed  independently  of  one  another.  Also  note  that  the  (potentially)  high-dimensional 
table  of  scores  need  not  be  stored  in  computer  memory,  since  the  maximum  score  can  be  computed  on  the 
fly,  thus  circumventing  any  memory  exhaustion. 

A  question  which  may  arise  while  implementing  the  proposed  algorithm,  is  whether  all  the 
polynomial  coefficients  should  have  the  same  discretization  resolution,  and  if  not,  how  to  determine  it. 
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From  Eq.  (2)  it  can  be  seen  that  the  small  change  A ak  in  the  kth  coefficient  (due  to  its  discretization 
resolution)  results  in  the  following  change  in  6  at  the  arc  length  s: 

A  0(s)  =  A  aksk  (A.l) 

Ideally,  we  would  like  a  uniform  resolution  for  0,  which  would  mean  that  A 0(s)  needs  to  be 
independent  of  s  and  k.  Although  this  dependency  cannot  be  eliminated,  it  can  be  minimized  by  choosing 
a  specific  value  for  A ak.  Assuming  the  desired  constant  value  S  for  A 0(s),  we  minimize  the  following 
squared  error  integral  to  obtain  the  optimum  value  A ak: 

r  Lyi XCLX 

A a*k  —  argmin  I  (|A0(s)|  —  S)2ds 

Aak  J-Lmax 

C  Lmax  2 

=  argmin  I  (|Aaksk|  —  S )  ds  (A.2) 

Aak  J-Lmax 

S  /  1  \ 

By  choosing  such  resolutions  for  the  polynomial  coefficients  ak,  and  similarly  for  bk ,  the  curve  space 
is  discretized  more  uniformly,  hence  increasing  the  accuracy  of  the  search  for  the  high  scoring  curves. 
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Fig  1.  (Left)  Different  possible  curves  passing  through  a  seed  point  are  tested  and  their  scores  are  computed. 
(Middle)  The  curve  with  the  highest  score  is  selected.  (Right)  The  process  is  repeated  for  all  the  remaining 
seed  points. 
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Fig  2.  Curves  starting  from  the  seed  point  x0  are  parameterized  by  the  arc  length,  s  E  [— L_,L+].  The  unit 
tangent  vector,  t(s),  is  approximated  with  polynomials. 


19 


•  *»'*$*  *»••*>. 
M«««n«UMn  l«i  * 

•  »*••«  ?  +  *  t  •  *  •  •••****»  st  •*€%*•»  e  *  «  * 

C  M  M  n#  |  •  5  i  D  M  O  3  •  M  M 

•  «  «  »  • 

•  *4$  ••  «  a 

«  *  •%•••«%  ♦  •  ♦  •  m  ft  »jf#*#*4«<»4»  Mr  #  m  • 
*#4#4*4#%4+«#4#*«**fc#9#4%«  pr*  *  •  »  •  • 

jt# 

•••••*#•+*  ♦  »*•+%•  s  #  ^  i«  *  <  »  +  •  ir  •  #  *  • 

•  »  s  •  •  *  »  *  1  •  «  #  #  *  £  ♦  I  ♦  %  »  ♦  »  •  4  «  •  t  9  4  % 

•  MMf#lMt«iliiitaaj#ilitia««a#» 
»•*%%••«•>  4#  *9  •’$'%«#•  .*■» 

f  m  n  •  #on  H  3  ’  ■<  ;  V,  ♦  r  nwH  •  *  !  •  •  m 

<  ♦  *  C  C  j  ♦  •  M  >  ?  >  M  H  i  •  •  4  m  ♦  •  *  •  «  t  •  •  <r 

•  4  •  •  %  *  ft  %  81  %  4  •  3  JP -*  %  4  ♦  ft  \  \  4  •  ®  ♦  •  4  ♦  #  *  ft  •  ft 

S!  i  •  }  M  <  «  M  •  -f  f  M  •  t  •  t.  «»  \  V»  •  M  3  •  ♦  t  ♦  »  • 

aaf i#l\i«aia«t#»aa 
99%  •  $  %  *  %  «i»«m  9  t 

444  •  **S*«$%4C«  *#*»#■» 

m««  2  233%  •  «  |a*«  !><•  •  » 

•  4  «  ft  4T  #  &»  4  ft  »  ♦>♦»♦! 

•a  •'  *4V  •#»•%•*  •  T  «k«.«0«>00» * 

•  . . 


Fig  3.  Reconstructed  ODFs  (top  left)  and  the  tractography  results  (rest  of  the  subfigures)  on  the  excised  rat 
spinal  cords,  using  various  values  for  the  bias  parameter  X. 
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Fig  4.  Tractography  results  on  a  human  brain  HARDI  dataset  using  polynomial  orders  of  (left)  N  =  3  and 
(right)  N  =  2,  shown  in  (top)  sagittal,  (middle)  coronal,  and  (bottom)  axial  views. 
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Fig  5.  Stereoscopic  rendering  of  Fig.  4  (bottom,  left).  To  see  this  image  in  3D,  please  cross  your  eyes  and 
move  the  image  closer  or  further  away  from  you  until  you  see  what  appears  to  be  a  third,  3D  image  in  the 
middle. 
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Fig  6.  Tractography  results  on  a  monkey  brain  HARDI  dataset  shown  in  axial  (top)  and  tilted  (bottom) 
views.  Seed  points  were  randomly  generated  inside  the  transparent  blue  regions  (top,  left).  A  mask  of  the 
tracts  is  shown  superimposed  on  the  T1  anatomy  image  (top  right  &  bottom). 


Fig  7.  Tractography  results  from  five  human  brain  HARDI  datasets  combined  using  geometric  (top  row)  and 
arithmetic  (row  two)  means,  and  from  individual  subjects  (two  bottom  rows),  shown  in  coronal  (left)  and 
sagittal  (right)  views.  24 


